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The Gaussian Effective Potential (GEP) is shown to be a useful variational tool for the study of the 
magnetic properties of strongly correlated electronic systems. The GEP is derived for a single band 
Hubbard model on a two-dimensional bi-partite square lattice in the strong coupling regime. At 
half-filling the antiferromagnetic order parameter emerges as the minimum of the effective potential 
with an accuracy which improves over RPA calculations and is very close to that achieved by Monte 
Carlo simulations. Extensions to other magnetic systems are discussed. 
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I. INTRODUCTION 

The Gaussian Effective Potential (GEP) is a well established variational tool Q, H, 0, 13, H, EH, 0, H, H which 
has been mainly used and developed for describing the breaking of symmetry in the standard model of electroweak 
interactions [13, [Tl|, H, E3, El, E3]- Quite recently, the GEP has also been shown to provide a very useful non- 
perturbative method for the study of condensed matter phenomena like superconductivity E3, and we expect 
that many other interesting phenomenological aspects could be described by the same tool. Quite generally, the 
method allows to study the broken-symmetry ground state of a bosonic theory describing quantum and thermal 
fluctuations for a large class of physical systems. Actually, a Bose field emerges in the theoretical description of many 
physical systems like superconductors, superfluids, magnetic materials, disordered fermions and many others. Thus 
we believe that the whole capabilities of the GEP have not been explored yet. 

In this paper we show that the GEP can be used for describing the strong-coupling limit of magnetic materials, 
where a spontaneous breaking of symmetry occurs at low temperature, and the local order parameter (i.e. the local 
magnetic moment) is a Bose field with excited states known as spin waves. In the strong-coupling limit, the usual 
perturbative methods are not reliable and the more important features are usually recovered by non-perturbative 
techniques like RPA[l9|. fluctuation exchange (FLEX) approximation |20l |. functional renormalization group (RG) 
[1H and numerical Monte Carlo calculations [12] ■ When there are no exact results available, Monte Carlo outputs 
are often regarded like experimental data even if such numerical calculations are plagued by several shortcomings 
like border effects, small size of the samples and a finite temperature. We show that the GEP provides an analytical 
variational tool for describing infinite systems at zero temperature with an accuracy which improves over RPA and is 
close to that achieved by Monte Carlo. 

In order to compare the results we study the well known half-filling Hubbard model in two spatial dimensions on 
a square lattice. The model has a broken-symmetry magnetic ground state at any coupling, and in the strong-coupling 
limit the local magnetic moment is known to saturate at a smaller value than predicted by the simple mean-field theory. 
According to Monte Carlo data [22] and RPA calculations the local moment is reduced by almost one half as 
a consequence of quantum fluctuations, in agreement with the theoretical limit which can be extracted from the 
Heisenberg model [HI, HH- The GEP predicts the correct strong-coupling limit of the Hubbard model, and we believe 
that such a variational method may provide a useful analytical interpolation for any coupling. Unfortunately we were 
not able to get the GEP exactly from the Hubbard model, and we had to approximate the exact spin wave Lagrangian 
by a power expansion. Neglecting higher order powers can be shown to be reasonable in the strong coupling limit, 
but there is no control of the approximation in the small-coupling regime and eventually at the quantum transition 
where fluctuations become very large. Thus we could not check the validity of the GEP in the small-coupling regime 
where, on the other hand, mean-field and perturbative methods are known to work well. That shortcoming is a limit 
of the present calculation while we expect that the exact GEP should give the correct prediction for any coupling. 
However the calculation has its merits and shows the potentiality of the GEP in the non-perturbative strong-coupling 
regime where alternative analytical methods are particularly welcome. Once the method has been tested on the well 
known Hubbard model, and its limits are understood, we believe that the GEP may represent a valid alternative tool 
for the study of other more complex magnetic systems. 

The paper is organized as follows: in the next section the GEP approximation for the Hubbard model is discussed; 
then the effective potential is evaluated at half-filling in section [Hi] where the phenomenological predictions of the 
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method are compared with the available Monte Carlo data. 



II. THE HUBBARD MODEL: DEFINITIONS AND GEP APPROXIMATION 



Let us define the single-band Hubbard model on a two-dimensional bipartite square lattice according to the Hamil- 
tonian 



H = - t zZ U{r f )c r7 (r>)+h.c.)+U , Y,n t (r)n l (r). 



(1) 



where cj T (r),c cr (r) are the usual fermionic creation and annihilation operators, (,) restricts the sum to nearest neigh- 
bour sites, and n a (r) = cl(r)c a (r) is the particle number operator at site r for spin a (Below we set U=U'/3 in order 
to agree with a widely used notation[l9l [22T|). 

According to the standard path-integral approach (26|. in the basis of coherent states for fermions, the fluctuation 
amplitude Z reads: 



D[W] 



jifdtL 



(2) 



where \& are the Grassmann variables associated with the fermionic annihilation and creation operators and L is 
the Lagrangian density. An auxiliary bose field 4> is introduced at any space-time point by the Hubbard-Stratonovich 
transformation, and the Grassmann fields can be integrated out exactly yielding 



with the effective action 



S eff(<f>) = - g 



Z = J D[(f>}e iS 'ff^\ 
dt ^2 <? *) - * m Dct ( id t +fJ,-M(<j>) 



written in terms of the matrix M(<f>) 



(f,t,a\M{ct>)\?',t',P) = 6 ttt , 



~2tS a ,j3 8?i>f + VU(f>(r,t) ■ Ta,f}Sr,r' 
(?") 



(3) 



(4) 



(5) 



where f a p are the Pauli matrices. The bosonic effective action S e ff is exact, and describes the spin- wave excited 
states of the Hubbard model at zero temperature[27| . 

We would like to build the GEP for the effective action Q , and as a first step we take a shift of the bosonic field 
and write = 0(r, t) as the sum of a background constant (non homogeneous) field 0o plus a fluctuating field <f) ' 



(f,t) = cf> (r f ) + 4>'(r,t). 



The effective action becomes 

Seff[4>0 + 4>'] 



j dt zZ \ + £') 2 - 1 lnDct { idt + ^ - M (^°) - SM (^')) 



where the matrix SM is defined according to 

(f, t, a\5M(<£')\r', t', (3) = 5 t ,t>VU(f' ■ t q , p 5?^. 
Then we can split the effective action as S e f / = So + 5int and write 



o 2 -HTrln UG i<; 



s * = -\j dt zZ 

r 

Sint = - J dt J2 -\j dt zZ 



(6) 
(7) 

(S) 
(9) 



Tr [iGo(<M8M{<t>') 



(10) 
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where iGo(4>o) is the mean-field fermion Green function in a background magnetic field 4>q: 

*G o (0o)= ( id t +A*-M(&)) ~\ (11) 

The expansion is exact provided that all the infinite terms are included. 

It is important to note that symbols here have a different meaning compared to the standard expansion reported 
in textbooks . Actually, in our calculation the background field <po is not fixed at the saddle-point of the effective 
action in ((H) , but it represents a variational parameter to be evaluated by a minimization of the Gaussian effective 
potential which includes quantum fluctuations. Similarly, the shifted field <p' does not describe just the fluctuations 
around the mean-field value of 4>o- it embodies the behaviour of quantum fluctuations for every different value of the 
background field. Thus in this variational calculation a different quantitative result is obtained by use of the basic 
standard formalism. 

Let us introduce the trial Gaussian action term Sgep 

Sgep = \Y,j dt4> , a (r,t)(g- 1 r b ^ b (r\t') (12) 

r,r' 

with the matrix elements g a b(r, 1, r', if) playing the role of variational parameters, and let us write the effective action 
as 

Seff = Sgep + (So + Smt — Sgep) ■ (13) 

According the amplitude Z reads 

Z[$q]= f D[4>']e iSGEP e l(Sa+s ^ t ~ SGEp] (14) 



and we can write 

lnZ$, ] = lnZ [&> ] + ]n(e^ So+s ^- SGE ^) GE p (15) 

where 



M'Po 



D[cj)'}e lSGEP (16) 



and with the Gaussian averages defined as 



fD[$ >}(... )e^ EP 
{{■■■))gep- — 7~^777\ - g ■ ( u ) 

plOGEP 



ID 

As usual [IcL . [l8|. the GEP is defined according to 

Vgep[(Po } = y (ilnZ [<f)o ] - (S + S int - S G ep)gep^ ■ (18) 

where V is a total space-time volume. The effective potential Vgep is the energy density of the system, and it 
obviously depends on the choice of the matrix elements g a b that will be regarded as variational parameters. The 
evaluation of the GEP then only requires the calculation of the gaussian averages in lfT8]) according to Wick's theorem 
with 

(cf>' a (r,t)(t>' b (r',t)) GE p=ig ab {f,t,r',t'). (19) 

Unfortunately the expansion in lfT0|) can only be evaluated up to some finite order, which means that some extra 
approximation is added besides the genuine variational GEP of ifl8|) . Neglecting higher order powers is a very 
reasonable approximation in the strong coupling limit, but is expected to fail in the small coupling regime. In fact the 
approximation may be regarded as an expansion in powers of the fluctuating field <\>' around the generic point (j>o- at 
any finite order, l|10p fits the exact effective action around 4>o in a neighbourhood that will get larger as more terms 
are added to the expansion. According to <fl8|) . the GEP is evaluated by a sampling of Si n t[<j>'\ through the Gaussian 
functional by (fT7|) . then the accuracy of the GEP depends on the width of the Gaussian functional e?q>(iSGEp)'- a 
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narrow trial functional will average the approximate effective action around <po on a small range where the expansion 
is a good fit of the exact effective action; a broader functional will average the effective action on a larger range where 
the expansion departs from the exact action, resulting in a totally uncontrolled approximation and wrong predictions. 
But the width of the wave functional is determined by the relative weight of fluctuations which, in turn, are related to 
the Coulomb coupling U ; actually, it is well known that large values of the coupling suppress fluctuations of the local 
moment (i.e. the Hubbard-Stratonovich bose field). In this sense, we expect that in the strong coupling limit even a 
second order expansion in l[l"0|) should be enough while in the weak coupling regime and at the quantum transition, 
where fluctuations become very large, the method is doomed to fail. However, in principle, the method could be 
improved order by order, even if the real calculation of the gap equation is expected to become very difficult. 

With such limitations understood, we can easily evaluate the GEP by use of the simple second order approximation, 
and content ourselves with the strong coupling regime of the model. Neglecting costant terms, we can write 



and the averages read 



hxZ [<fo] = ^Trln<? (20) 



(S )gep = -\ J dt^^ir) + iTrln^Go^o)) (21) 



(S m t) GE p « -\ J dtY,($'\r,t)) GEP - ^Ti-((iGo(M6M($')) 2 )gep (22) 
f 

where we have summed up to second order and have omitted odd terms whose average is zero. We notice that 
the variational parameters g ab are implicit functionals of the local background field <fio(r), as they are going to be 
determined for any local configuration of the field by requiring that the effective potential is stationary. Using ([8j 
and (fl9l) we obtain 

T±{(iG {fo)$M{$ l )) 2 )G E p = iU JdtJ dt'J2(G 0aa ,(r,t;r',t'J )T«, fJ x 

r,r' 

x Gop/3'(r",t';r,t,(j>o)T^ a g ab (f',t';r,t)). (23) 
Finally we can Fourier transform and according to (fl8j) we are left with the effective potential 

Vgep [<M = \ V / In g aa (q, 0) - i V / In (iG 0ota (q, fi, $ )) + i- f dtJ2 $o 2 (r)+ 

+ \ E / 9aa(Q, A) - | / 9ab(q, tt)A^(<f, fi) (24) 

where and 



K ab (q,Q) = Tr / G {k,u, Mr a G (k - q,u - n, <j> Q )r b . (25) 

For any local configuration of the background field 4>o(r) the matrix g follows from the stationary condition 
SVgep I '5g a b = (gap equation) which reads 

(g- 1 )^ = -5 ab -iUK ab . (26) 



This gap equation defines the matrix elements g a b as implicit functionals of the field tfio (f) , and by insertion in ([18] 
the GEP becomes 



Vgep[M = / V9 aQ (g,tt)-*V I ln(iG 0aa (q,n,$ o j)+^- f dtV0 o 2 (r). 



(27) 



where we have omitted all constant terms. The three terms in (|27|) are easily recognized as the quantum spin 
fluctuation energy, the electron energy and the classical spin energy respectively. 
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A closer inspection of (|26j) allows us to obtain a deeper understanding about the nature of our approximation and 
its relationship with other approaches to the Hubbard model. From a formal point of view, (|26|) tells us that the 
spin correlation matrix g can be seen as the sum of a geometric series whose common ratio is iuK. This leads to 
a diagrammatic expansion of g in terms of bubble and ladder graphs typically involved in the RPA calculation of 
quantum fluctuations. So, the GEP does not introduce the resummation of any different class of diagrams as, e.g., the 
FLEX approximation [20j]. Anyway, RPA and GEP are not equivalent. In fact, in the random phase approximation 
the evaluation of quantum corrections to the sublattice magnetization makes use of the one-particle Green's function 
obtained from the mean-field value of the background field 4>o Oil- In the GEP, instead, as already pointed out 
above, the background field dependence of the fluctuation matrix relies on a variational minimization, so avoiding to 
fix </>o at its mean-field value. Furthermore, the reliability of our approach in the str ong and very strong coupling 
limit makes it a complementary tool to RG methods like the Interaction Flow scheme [21(, which become inaccurate 
when the interaction is equal or greater than the bandwidth, but works very well for small values of U/t. 

The ground state magnetization is determined by the stationary condition of the effective potential 6Vgep/5<I>o = 0. 
The mean-field result is recovered by neglecting the quantum fluctuations: in fact if we neglect the first term on the 
right hand side of l|27p the stationary condition reads 



00 (f) 



UTi 



tG 



o(<Po 



(28) 



which is the standard mean-field self-consistency equation. The functional dependence of the matrix elements g a b on 
4>o ensures that the minimum of the total effective potential occurs at a shifted value with respect to the mean-field 
magnetization. 

Furthermore, (|28| establishes a relation between the fluctuation field and the sublattice magnetization. In fact, 
since the local mean-field magnetic moment is equal to 



(S(r,t))=Tr 



-G (</>o) 



(29) 



where the spin operator is 



S(r) = 4(r)^f-c CT (r), 



(30) 



we can easily write 



(31) 



This relation is more general than mean-field approximation and can be also recovered from the definition of the field 
$ (|H|)- Keeping in mind this result, we will consider the effective potential Vqep as a function of the magnetization 

m=2\(S(r,t))\. 



III. RESULTS AT HALF FILLING 

The half-filled Hubbard model is known to be an antiferromagnet for any strength of the coupling U. We assume 
a staggered magnetization on the bipartite square lattice 



(j>o(r) = ±ipn 



(32) 



where n is any constant unimodular versor, |n| = 1, and the sign is positive on a sub-lattice and negative on the other 
one. The effective potential becomes a function of the scalar ip which gives the strength of the local magnetization. 
The versor n breaks the rotational symmetry of the model, and the opposite signs of 4>o on the sub-lattices break the 
translational symmetry of the square lattice. The unit cell doubles its size as it is replaced by the unit cell of the 
sub-lattice. The first Brillouin zone reduces to one half and its boundary is at the Fermi surface of the unperturbed 
electron gas at half filling (perfect nesting) . Let us define the two-component Green function G ± (ip) 



uj±e{k) 



1 a /3 



lu 2 — E 2 (k) + irj 



(33) 



where the free-electron band energy e(k) is 



e{k) = -2<^]cos(fc a ) 



(34) 



and the mean field band reads 



E(k) = Ve 2 (fc) + £V 



(35) 
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Figure 1: Effective potentials in units of the bandwidth t for y=10 (a) and y=30 (b). The minima of the effective 
potentials, i.e. the ground state magnetization, occur for m=0.55 in (a) and m=0.6 in (b) (see further details in the text). 



At half filling (/i = 0) the electron Green function Go ((f) reads|26] 

Go(k,u;ip) = G + (k,uj;ip) 



(36) 



for any k belonging to the reduced first zone (below the unperturbed Fermi energy), while outside the reduced first 
zone (above the unperturbed Fermi energy) the Green function is given by the second component 



G (k,uj;ip) = G (k,w;<p). 



(37) 



Insertion into f25j) and (|26f yields the spin wave correlation matrix g. 

The ui integration in lj25"|) can be carried out analytically. In this way we are left with a two-dimensional k-integration 
for the real part of g and a unidimensional one for its imaginary part, both to be performed by numerical methods. 

The real part has been obtained integrating by Simpson's rule on a 70 x 70 grid in the full first Brillouin zone, while 
for the imaginary part we used a 3000 point sampling of the interval [—it, 71"]. 

The effective potential Vgep can now be evaluated through lj27|) . Exploiting the parity properties of g(q, fi) the 
first integral in lj27|) can be restricted to positive values of Q, and to the upper right quarter of the Brillouin zone. The 
employed grid contains 200 x 20 x 20 points in the (O X q) space. 



Figu re Q] p resents two plots of the effective potential in units of the bandwidth t for y =10 (figure [1(a) I and y=30 
(figure 1(b) ). As expected, the effective potential shows a minimum for a non vanishing value of to, which corresponds 
to the magnetization of the broken symmetry ground state; a more accurate estimate of to is then obtained by a least 
square fit of the region of the minimum with a parabolic curve. 

Another important feature of our approximation can be observed in figure[TJ the potential broadens as -jr decreases. 
This behaviour implies that the smaller is the coupling the larger are the fluctuations of the field, so we expect that 
the simple second order expansion of lfT0|) should fail in the small coupling regime. 
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Figure 2: Sublattice magnetization mvs j within the GEP (blank squares), RPA (stars), Monte Carlo (errorbars) and 
mean-field (crosses) approximations. The arrow indicates the spin-wave strong coupling limit m=0.606. 



This trend is evident when we compare our results for the magnetization with other calculations present in literature. 
As shown in Fig[2j for y<4, the RPA [iji] and Monte Carlo [22] predictions converge reasonably well towards the 
mean field approximation, while the GEP gives wrong results due to a bad averaging process by a Gaussian functional 
which becomes too wide, as discussed in the previous section. However the approximation improves in the strong 
coupling region (~ ^8). Data show a tendency to saturate to a limit value and this trend appears slower than with 
RPA, in agreement with Monte Carlo results. In order to get the satura tion v alue we performed a calculation using the 
very strong coupling constant y=30, obtaining m=0.6 (see also figure 1(b) |. This limit represents a very important 
consistency test for the GEP. It is in fact straightforward to show that for large j the Hubbard Hamiltonian (fl]) reduces 

to the anti-ferromagnetic Heisenberg model with exch ang e coupling A comparison between our magnetization 
and Monte Carlo simulations for the Heisenberg model [25| points out that the GEP has the correct limiting behaviour 
(figure [3]) . This result is also in agreement with spin- wave theory predictions [HI- Moreover, also in the intermediate 
region our results nicely interpolate between the weak and the strong coupling limits and are closer to Monte Carlo 
calculation than RPA. 

In conclusion, we have shown that the GEP can be regarded as a useful non-perturbative tool for studying magnetic 
systems with a broken symmetry ground state, and that it can be regarded as a complementary tool together with 
other well established approaches like functional RG, FLEX or RPA. In particular, we have been able to evaluate 
the effective potential for the Hubbard model and to write down an analytic expression ( ([27]) ) where the "classical" 
electron and spin energies and the quantum spin fluctuations are clearly recognizable. 

For the half filled Hubbard model we compared the magnetization with the result of other approximations showing 
that the GEP provides an improvement over the RPA approach in the intermediate and strong coupling limit, while it 
seems to be inaccurate for small values of — . However in the weak coupling limit the failure seems to be a consequence 
of the simple second order expansion, and we expect that the exact GEP should predict the correct magnetization. 
In this limit an improvement could be achieved by insertion of higher order terms in the expansion IjlOp . 

Further details about the GEP approximation can surely be obtained investigating its dependence on the band 
filling. The method of section Hill still holds for any filling of the band, while the half-filling electron Green's function 
(|33"ll should eventually be replaced by the general one which depends on the filling of the band. That would allow us 
to construct a phase diagram for the 2D Hubbard Model and discuss the competition between antiferrommagnetic 
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Figure 3: magnetization limiting behaviour for the GEP (full circles) as a function of the inverse coupling constant (U/t) -1 . 
The straight horizontal line shows the Heisenberg Model magnetization as obtained with a Monte Carlo (MC) by Sandvik [25l ] 
while the triangle represents the spin-wave result for the same model [24l |: they are shown in correspondance to (U/t)^ 1 as in 
that limit Hubbard and Heisenberg model should coincide. As it is evident from the figure, the GEP magnetization trend 
confirms this expectation, supporting our claim for the GEP to be very reliable in the strong coupling limit. 



order and d-wave superconductivity[13, HH 
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